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I. INTRODUCTION 

Topological insulators represent a new class of mate- 
rials where the topology of the bulk electronic structure 
induces non-trivial effects such as the emergence of edge 
states. The edge states are robust against smooth de- 
formations of the crystals or the presence of disorder. ^'^ 

Among all known classes of topological insulators, the 
time-reversal invariant ones have a special status be- 
cause they were already engineered and characterized in 
laboratories.^"^ As pointed out early in the development 
of the field, the bulk-edge correspondence in time-reversal 
invariant insulators obeys a Z2 classification.^ It was also 
clear from early stages that the topology of the time- 
reversal invariant bulk electronic structures is classified 
by the twisted Real K-Theory,^ which pointed again to a 
Z2 classification. Homotopy arguments lead to the same 
conclusion. 

Several explicit formulations of the Z2 invariants were 
given along the years but all of them involve globally 
smooth gauges. Gauge independent invariants were for- 
mulated in Ref. 11, but their effectiveness remains to 
be tested in 3 dimensions. For time-reversal invariant 
systems, which inherently have trivial Chern numbers, 
an important result by Panati assures the existence of 
such globally smooth gauges. Unfortunately, the result 
by Panati is not constructive and at this point we don't 
have a standard algorithm to construct globally smooth 
gauges, something that in many instances proved to be 
a formidable task. 

The early work of Ref. 3 proposed to look at the Pfaf- 
fian of a particular overlap matrix as function of the k- 
vector. Generically, this Pfaffian can become zero at iso- 
lated /c-points, which always come in pair. The number 
of paired first order zeroes, taken modulo 2, was found 
to be a topological invariant. Several equivalent expres- 
sions of the Z2 invariant were derived in Ref. 13. This 
work introduced the notion of time-reversal polarization, 
which was shown to be quantized modulo 2. The compu- 
tation of the time-reversal polarization requires a globally 
smooth gauge, which must also be adapted to the time- 
reversal symmetry (see Eq. 3.10). The time- reversal po- 
larization approach inspired yet another formulation of 
the invariant, involving the Pfaffian and the square root 



of the determinant of another overlap matrix, computed 
at the time-reversal invariant /c-points. This formula- 
tion played a special role since it admitted extentions to 
higher dimensions. ^'^^ Furthermore, the Z2 invariant was 
formulated as an obstruction against achieving a glob- 
ally smooth gauge of certain kind, leading to yet another 
equivalent expression involving the Berry curvature and 
the Berry phase. Starting from this later expression of 
the invariant, Fukui et al were able to develop a (almost) 
gauge-independent method of calculus in 2-dimensions 
(the method still requires a time-reversal adapted gauge 
at the boundary of half of the Brillouin zone).^^ This 
method was later extended in 3-dimensions,^^ and it be- 
came the method of choice when computing the Z2 in- 
variants for non centro-symmetric systems. 

The requirement of special smooth gauges in the clas- 
sic formulations of the Z2 invariants is unfortunate, and 
the (almost) gauge-independent method of calculus de- 
veloped by Fukui et al can be quite involved. For exam- 
ple, the prediction of the first topological insulators in 
3D was in great part possible because the Z2 invariants 
simplify tremendously when inversion symmetry is also 
present.^ Without inversion symmetry, the evaluation of 
the Z2 invariants remains a very difficult task. For ex- 
ample, in a study on strained bulk HgTe material, a 
non centro-symmetric system, even with a tight-binding 
model it was more convenient to complete slab calcula- 
tions and look directly at the surface states rather than 
compute the bulk Z2 strong invariant. A similar ap- 
proach was followed in a recent first-principle study on 
the non centro-symmetric metacinnabar compound. So 
evaluating the Z2 invariant is already difficult at the 
level of tight-binding modeling, but the difficulty be- 
comes overwhelming when attempting first principle elec- 
tronic structure calculations. This aspect was recently 
discussed in Ref. 25, where a solution was proposed based 
on hybrid Wannier functions. Subsequent work,^^ has 
also employed hybrid Wannier functions to derive equiv- 
alent Z2 invariants in 2 dimensions. Notably, this later 
work made use, like us, of the full (not just the trace) 
adiabatic connection. The use of hybrid Wannier func- 
tions to efficiently re-formulate the topological invariants 
was originally introduced in Ref. 27. 

In 3 dimensions, topological insulators were also shown 
to display quantized magneto-electric polarization,^^ 
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which can be written as a Chern-Simons integral. This 
invariant was shown to be completely equivalent to the 
previously introduced strong invariant. The com- 
putation of the Chern-Simons integral requires again a 
globally smooth gauge and the fundamental difficulties 
introduced by this requirement were already highlighted 
in Ref. 30. To date, nobody has achieved a direct evalua- 
tion of this Chern-Simons integral, even for tight-binding 
models. The magneto-electric polarization was computed 
indirectly, using the second Chern number and dimension 
reduction technique. ^^'^^ It became a sure fact that com- 
puting the second Chern number in 4 dimension is much 
easier than computing the Chern-Simons integral in 3 di- 
mensions, and this is precisely because the second Chern 
number admits a manifestly gauge indepent expression 
based entirely on the projector onto the occupied states. 

Our present work provides equivalent formulations of 
the Z2 bulk invariants that are manifestly gauge inde- 
pendent. We report results for both 2 and 3 dimensions. 
The new formulas are computationally trivial for both 
tight-binding and first principle approaches. The key to 
these results is a "monodromy" argument, which was re- 
cently introduced in Ref. 32. Basically, instead of looking 
at the polarization, we examine the full non-abelian adi- 
abatic transport along time-reversal invariant lines in the 
Brillouin zone and take advantage of the special behavior 
under the time-reversal operation. Using the elementary 
properties of the determinants and Pfaffians, we are able 
to show that the determinant of the monodromy, com- 
puted along closed time-reversal invariant paths in the 
Brillouin torus, can be written as the square of a well de- 
fined quantity. This quantity divided by the square root 
of the determinant of the monodromy takes the quan- 
tized values of ztl, and becomes the building block for 
our invariants. 

In 2 and 3-dimensions, we look at pairs of time-reversal 
paths on the Brillouin torus. For such pairs, we show that 
the square root of the determinants of the monodromies 
can be taken in a canonical way, allowing us to define 
a true topological invariant for each such pair. In 
2- dimensions, this construction gives the unique Z^ in- 
variant, while in 3-dimensions it generates four indepen- 
dent weak invariants, which can be subsequently used to 
generate the strong Z^ invariant. 

We use tight-binding models with time-reversal sym- 
metry to test our formulations and to show how the con- 
struction works. The models include interactions that 
strongly break the inversion symmetry. 



II. THE MAIN CONSTRUCTION 

Let us consider a one dimensional, translational in- 
variant lattice model, described by a Bloch Hamiltonian 
iJ/e, di N X N complex matrix with /c-dependent entries. 
The Hamiltonian acts on the fixed space of TV-component 
complex vectors, the space. We assume time-reversal 
symmetry, that is, we assume the existence of an antilin- 



ear operator acting on C^, such that: 

OHkO-^=H_k- (1) 

We also assume that we are dealing with an insulator, 
so Hj^ is assumed to have a spectral gap at all /c's. We 
denote by Pk the projector onto the states below this 
spectral gap. 

Our construction starts from the monodromy Uk,k' de- 
fined as the unique solution to the following differential 
equation: 

i£Uk,k' =i[Pk.dkPk]Uk,k', (2) 

with the initial condition Uk' ,k'=Pk' - Here, k' is an ar- 
bitrary but fixed /c-point. The monodromy provides a 
parallel transport, that is, an isometric mapping of the 
space Pk'C^ into the space PkC^ , more precisely: 

Pk = Uk,k'Pk'U-l. (3) 

The monodromy is also known to generate a one param- 
eter unitary group: 

Uk.k'Uk'^k" = Uk^k", Uk^k'Uk'.k = Id. (4) 

Eq. 2 can be derived, and it was first derived (see 
Ref. 33), from the modern formulation of the Adiabatic 
Theorem. If a local gauge (i.e. smoothly k-dependent 
bases for PkC^ spaces) was pre-defined, then Eq. 2 takes 
the more familiar form: 

^U{k) = iA{k)U{k), (5) 

where A{k) is the full non-abelian adiabatic connection 
discussed by Wilczek and Zee in Ref. 35. For a more 
detailed discussion one can consult Ref. 36. We will, 
however, want to stay way from the later Eq. 5 because a 
smooth gauge can be, in general, chosen only locally. And 
even though time-reversal invariant systems admit global 
smooth gauges, constructing such a globally smooth 
gauge can be quite a formidable task. 

One key observation is that Eq. 2 can be integrated 
without making use of any gauge. Indeed, assume that 
we divided the interval [A:', k] in small subintervals: k' = 
ki . . . < kn = k^ and let us form the product: 

Ui = Pk,Pk,_,...Pk,- (6) 

By a simple term counting, one can easily see that Ui 
satisfies the equation: 

(7) 

= i{{Pk, - Pk,_,)Pk,_, - PkXPk, - Pk,-^)}Ui-i 

But this equation is nothing else but the finite difference 
version of our original Eq. 2. In other words, Eq. 2 can 
be integrated by forming the sequenced product shown 
in Eq. 6, using a fine-enough finite difference step. This 
discussion is not limited to one dimension but it can be 
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applied to the parallel transport along any arbitrary path 
in higher dimensional Brillouin zones. On a more tech- 
nical note, let us state that the projectors can be 
computed without using any particular gauge. There is 
quite a substantial number of different ways to accom- 
plish that, but just for the sake of explicit ness, let us 
mention that the projector onto a particular eigenvalue 
ei{k) can be computed as: 



(8) 



where Fi is the so called interpolating polynomial defined 
hj Fi{ej) = 5ij. 

Now let us conjugate Eq. 2 by ^: 

^{^l^M'}^"' = 0{i[Pk,dkPk]Uk^k'}0-\ (9) 
which leads to: 

ij^{OUu^k'0-^)=i[P-k.duP-k]{OUu^k'0-^). (10) 

Also, OUk',k'0~^ becomes the identity on P-k'C^ . This 
leads us to conclude that: 



-k'- 



(11) 



Therefore, if we want to compute the monodromy from 
-TT to TT, we can write: 



^7r,-7r — ^7r,0^0,-7r 



(12) 



The monodromy t/7r,-7r maps the space P-kC^ into itself. 
We can therefore enquire about the determinant of this 
monodromy. It is a fact that the determinant of Ut^^-tt 
is always equal to one for systems with time-reversal and 
inversion symmetries. If the inversion symmetry is bro- 
ken, this determinant can take, in principle, any value 
on the unit circle. As we shall see in the following, for 
the most general situation, the determinant can be writ- 
ten as the square of a well defined quantity, which will 
become the building block for our Z2 invariants. 

To see this, let us choose an arbitrary basis in P^C^ 
and PqC^ , which we denote by {e^} and {e^}, respec- 
tively. Before we start the calculation, let us point the 
following fact (Q = arbitrary linear operator acting on 
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We can put the above fact in a more convenient form. 



(e°|^Q|e^) = E.(e°|^|eO)(eO|g|eO), (14) 

in which case we see a simple rule, that when inserting 
an identity operator J2s ^fter the anti-linear op- 

erator 0, all the resulting matrix elements after must 
be complex conjugated. 



We can now start the calculation: 

(e-|C7,,_,|e^) = {el\U^,oOU-le-'\e-p) 



(eS|t^.,o|eg)(eg|^|eO)(eO|C7-J|e|)(e||^-i|e^). 



(15) 



Summation over repeating indices was assumed above. 
We denote by U the matrix of elements 



Note that 



Uap = {el\U^,o\es)- 



(e°|[/-J|ep = [/-/. 



(16) 



(17) 



Also, since C/t^^q is an isometry, the matrix U is unitary 
and consequently: 



(eO|[7-J|ep = U^^. (18) 

We denote by Oq and 0^^ the matrices of elements: 

(OoU = (eyie'e), {O.U = {eZ\0\e],) (19) 

and we point out the identity: 

{el\0-^\e-p) = {e-')^p, (20) 

where the complex conjugation is due to the property 
stated in Eq. 14. With these technicalities behind us, we 
can now state that: 



Therefore: 



U^,-^ = UOoU^O-\ (21) 



det{[/^,_^} = det{UOoU^O-^} (22) 



and, since the matrices are antisymmetric, we can use 
their Pfaffians and the elementary properties of determi- 
nants to conclude: 



det{C7^,_^} = [Pf{i}-Met{C/}Pf{^o}]'. 
We arrived at our main conclusion: 



Pf{^^}-Met{£/}Pf{^o} 
^det{C/^ 



= ±1. 



(23) 



(24) 



The left hand side of Eq. 24 depends on the branch of the 
square root we chose, but once this choice is made, the 
value of the left hand side cannot be changed by smooth 
deformations of the Hamiltonian that keep the insulating 
gap opened. 

One can verify explicitly that the formula is completely 
independent of the bases chosen at k=0 and k=7r. In- 
deed, if we make a change of bases: 



(25) 
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then 



and 



det{;7} ^ det{W„} det{U} det{W^o}" 

Pf{4} ^ det{W,}Pf{4}, 
Pf{^o} ^ det{Wo}Pf{^o}, 



(26) 



(27) 



SO the invariance follows automatically. In fact, Eq. 24 
can be evaluated without making reference to any basis 
set, by just using the abstract (fundamental) definition 
of the determinant and Pfaffian.^^ 

Because we don't have a canonical way to choose the 
branch of the square root at the denominator in Eq. 24, 
we cannot assign a true topological meaning to this for- 
mula. For instance, it will be impossible to compare two 
separate systems, unless we have an explicit way to de- 
form them into each other without closing the direct en- 
ergy gap. This shortcoming can be eliminated in 2 and 3 
dimensions, where true topological Z2 invariants can be 
defined. This is discussed in the following sections. 

We should point out that, by using Eq. 23 inside the 
square root, Eq. 24 can be also written as: 



Pf{^^o} ( Pf{4} 



= ±1, 



(28) 



which shows the direct connection between our formu- 
lation and Eq. 3.24 of Ref. 13. The only difference is 
that we provide a different, but equivalent criterion for 
choosing the branches of the square roots. 



III. THE Z2 TOPOLOGICAL INVARIANT IN 
2-DIMENSIONS 

In 2-dimensions, we can follow a pair of time-reversal 
invariant paths on the Brillouin torus, and generate a 
pair of quantized numbers like in Eq. 24, such as: 



^0 



for the path 



and 



for the path 



_ Pf{^^(o,.)}-Met{?7o}Pf{^^(o,o)} 

\/det{t/(0,7r),(0,-7r)} 

k = (0, -tt) -^k = (0,7r), 
. Pf{^(.,.)}~^det{?7^}Pf{^V,o)} 

\/det{t/(7r,7r),(7r,-7r)} 



k = (tt, — tt) ^ = (tt, tt). 
We now form the product 

^2D = ^O^TT, 



(29) 

(30) 
(31) 

(32) 
(33) 



in which case the arbitrariness in choosing the branch 
of the square root at the denominators becomes irrele- 
vant because now we have a canonical way to chose the 
same branch for the square roots of det{t/(o,-7r),(o,7r)} 
and det{[/(7^ 7^)}. Indeed, the paths described in 

Eqs. 30 and 32 can be deformed into each other with- 
out breaking the loops or leaving the Brillouing torus. 
Therefore, the Bloch Hamiltonians i?(0, ky) and i^(7r, ky) 
can be adiabatically connected without closing the en- 
ergy gap, which means t/(o,-7r),(o,7r) can be continuously 
evolved into /7(7^ 7^), and same can be said for their 

corresponding determinants. Therefore we have an ef- 
fective way to make sure we choose the same branch of 
the square root for both determinants. It is totally ir- 
relevant which branch we chose (as long is the same) 
because if we change the branch for both square roots, 
then a minus sign appears twice and nothing changes. 
The conclusion is that E2D can be given a meaningful 
topological content, and different time-reversal invariant 
systems can be classified according to the corresponding 
value of E2D- The trivial insulator is contained in the 
class with = +1 and the topologically non-trivial 
insulators are contained in the class with ^20 = — 1- We 
could have started the entire construction from paths ori- 
ented along the kx direction, but this would have led to 
the same topological invariant (see the argument by Roy 
in Ref. 38). 

Let us follow right away with a non-trivial exam- 
ple. We chose to work with the Berne vig-Huges- Zhang 
model, ^ including the 5'z-nonconserving term discussed 
in Ref. 39. The model is described by the Bloch Hamil- 
tonians acting on the space: 



/ h{k) T{k) 
~ \ r{k)^ h*{-k) 



(34) 



where h{k)=d{k)-cr^ with cr={ax^(Jy^(Jz) and: 

d = (^dsin/ca^, Asin/cy, A — 25(2 — cos/ca; — cos/cy)). (35) 
The r term is given by:^^ 

sin kr — i sin k,, 



r(A:) = iK 







sm kx -\- 1 sm ky ' ^ ^ 



As written above, the model is symmetric to time- 
reversal and to inversion symmetry operations. We in- 
clude an additional term which will specifically break the 
inversion symmetry but leaves the time-reversal symme- 
try intact. To be as explicit as possible, let us mention 
that the action of the time-reversal operation = e^^^yR 
{K= complex conjugation) in is: 










1 


1) 

















II 




-1 

















-1 





I) 





(37) 
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The inversion operation is implemented by: 



P = 



(38) 



The additional term to the Hamiltonian that we consider 
here is: 



R 






Q _^iki-\-2ik2 
-ik\—2ik2 Q 





-ik\ — 2i/c2 












y g-*Aci-:.*/C2 / 

(39) 

The factor 2 in front of k2 was chosen just to introduce 
an anisotropy. R is the coupling constant. 

li R = 0, Hj^ displays topological phases for < 
A/5 < 4 and 4<A/5<8, and trivial phases for A/B<0 
or A/B>S.^^'^^ The insulating gap closes when A/5=0, 
4 and 8. For a finite R, the phase diagram changes; the 
energy gap closes at 4 points and a new topologically 
trivial phase appears. Let us be explicit and fix some pa- 
rameters, from now on, as follows: A = B = 1, A = 0.5 
and R = 1. Upon varying the parameter A, we found 
that the energy gap closes at 1.15, 3.34, 4.65 and 6.85. 
By just taking into account the known phase diagram at 
R = 0,^^'^^ it is naturally to assume that the topological 
phases occur when A is in between 1.15 and 3.34, and in 
between 4.65 and 6.85. Topologically trivial phases are 
expected in rest (see Fig. If). But this we are going to 
check explicitly. 

The first part of the calculation relates to the square 
root of the determinants of the monodromies. For this 
we have considered paths along the ky direction: 



k = {kx, -tt) k = (/Ca;, tt). 



(40) 



for which we computed the monodromy t^(/c^,-7r),(/ca.,7r)5 
by straight implementation of the Eq. 6, using 1000 dis- 
cretization points. We will be more specific about this 
part of the calculation shortly. We then plotted in the 
complex plane the value of the determinant of the mon- 
odromies as function of kx^ when kx varied from to tt 
(here we used again 1000 discretization points). Fig. 1 
shows the plots for diff'erent values of A. When taking 
the square root of the determinant, what we must have in 
mind is the Riemann surface of the complex function ^/z, 
shown in Fig. 2. Most of the available softwares, when 
given a complex number z in the plane, it will automati- 
cally place z on the top sheet of the Riemann surface. As 
explained above, we can chose any branch of the square 
root for the determinant at kx = 0^ but after that we 
must be consistent with this choice when we compute 
the square root of the determinant ai kx = tt. So we 
will always place the determinant at kx = on the top 
sheet of the Riemann surface. Then, by following the 
evolution of the determinant of the monodromy as kx is 
varied from to tt, we will be able to tell exactly where 



this determinant is located on the Riemann surface. If 
the determinant ends up on the top sheet, we don't need 
any correction, but if it ends up on the lower sheet, we 
must correct the output from the software by multiply- 
ing the square root by a correction factor a = —1. In 
Fig. 2 we chose several situations and explain in detail 
how a works. To summarize, in the actual calculation we 
let the software (in this case MATLAB) to compute the 
square root of the determinants and afterwards we cor- 
rected the square root of the determinant at kx = tt by 
the sign factor ce, which was determined from the inspec- 
tion of the graphs in Fig. 1 and using the prescription 
given in Fig. 2. While here we chose to use a visual 
inspection to determine a, and this was mainly to show 
the reader how things work, it is important to notice that 
a can be determined in an automated fashion, without 
any visual inspection. This observation is important for 
the implementation of the formalism in the first principle 
codes. 

The second step of the calculation consist of evaluating 
the Pfafiians and the monodromies at the nominators in 
Eqs. 29 and 31. These were numerically evaluated in the 
following way. The line between ky=0 and ky =7r (this 
time assuming kx to be either or tt) was discretized us- 
ing 1000 points and was diagonalized at all these ky- 
points. This provided us with four eigenvectors il^i{ky)^ 
i = 1, . . . , 4, sorted according to their eigenvalues. Being 
part of the space, the eigenvectors are represented as 
4-component column matrices. Note that the diagonal- 
ization procedure gives random phases for the eigenvec- 
tors, but this is irrelevant when we form the projector 
onto the first two egivenvectors: 

Pk, = \Mky)){Mky)\ + \Mkv)){^2(ky)\. (41) 

The projectors were represented as 4x4 matrices. As 
ky was progressing from to tt, we have continu- 
ously updated the monodromy matrix: U^PkyU^ start- 
ing initially from U=Pky=o- After the monodromy 
was computed, we used the bases {'0i(O), ^^2(0)} and 
{'0i(7r), '02(7r)} for the occupied spaces at A:^ = and 
TT to compute the 2x2 matrix U. 

The Pfafiians at A:^ = and tt are simply equal to 
(V^i(0)|6>|V^2(0)) and (V^i(7r)|6>|^2(7r)), respectively, and 
these matrix elements were easily computed using the 
action shown in Eq. 37. One important note here is that 
the Pfafiians have to be computed using the same bases 
(including the phases!) as the ones used to compute U. 
One can verified how converged the monodromies are, by 
examining the absolute values of their determinants. In 
general, these values will be less than 1, but converge to- 
wards this ideal value of 1 as more points are added to 
the discretization of the paths. 

We now return to Fig. 1 and discuss the results. We 
picked three A values in each region of the phase diagram, 
so that we have values close to the end points where the 
gap closes and values far away from these points. The ac- 
tual values are shown in the middle of each panel. Each 
panel shows, in the complex plane, the calculated value 
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FIG. 1. Results for the 2-dimensional model of Eqs. 34 and 39, with the parameters fixed at: A = B — 1, A = 0.5 and R — 2. 
Each panel shows the path in the complex plane of det{U(^kx,-7T),(kx,'^)} is varied from to tt, for different values of A. 

The panels are grouped into bundles of 3 (for example dl, d2 and d3), and each such bundle samples a region of the phase 
diagram where the energy gap stays open. In each panel, one can read the value of A, the correcting factor a and the value of 
the Z2 invariant Panel (f) shows the predicted phase diagram of the 2-dimensional model. 



for det{[/(/j.^ (/^^ 7^)}, as kx was varied from to tt. 
This determinant is always on the unit circle, but just 
for a better representation, we have artificially shifted its 
value inside the unit circle (by multiplying with the func- 
tion e ^), so that we can follow its intricate behavior. 
Given these curves, and assuming that the determinant 
at kx 

was on the upper sheet of the Riemann surface 
of the square root, we can easily determine the position 
of the determinant ai kx = tt on the Riemann surface of 
V^, and therefore the value of a (see the discussion in 
Fig. 2). We have placed these values directly inside the 
panels, together with the value of the Z2 invariant. Be- 
sides these calculations, we have performed calculations 
with a much more refined sampling of A, confirming the 
phase diagram shown in in panel (f) of Fig. 1. 



IV. THE Z2 TOPOLOGICAL INVARIANTS IN 
3-DIMENSIONS 

In three dimensions, we can use different pairs of time- 
reversal invariant paths and construct weak Z2 invariants 
first. Let us consider the following explicit pairs: 



pair 1 



{ 



k (0,0, -tt) k = (0,0, tt) 

k (0, TT, -tt) ^ /c = (0, TT, tt) 



and 



pair 2 



■{ 



(7r,0, -tt) -> fe = 
= (tt, tt, —tt) k 



(7r,0,7r) 
= (7r,7r,7r), 



(42) 



(43) 



for which we construct the corresponding 2-dimensional 
(weak) Z2 invariants, and ^2^, as described in the 
previous section. All we have to check, and this is ob- 
vious, is that the paths in each pairs can be deformed 
into each other continuously without leaving the Bril- 
louin torus (in fact, all four paths can be deformed into 
each other). The strong invariant is given by their prod- 
uct: 



EsD — '^2d'^2D' 



(44) 



We can start the construction from different pairs of 
paths, but at the end we can generate at most 3 indepen- 
dent weak invariants plus the unique strong invariant, a 
fact that can be shown by using a fairly general method 
introduced by Roy.^^ 

Let us again follow with an example. We chose to 
work with the model Hamiltonians reported in Ref. 41. 
We will tune the parameters for Bi2Se3, following Ref. 42 
(see Eq. 31 and Table II). Explicitly, we considered the 
Bloch Hamiltonians: 



Hk = 




(45) 
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FIG. 2. (Color online) The complex functionA/^ is multival- 
ued and its proper representation is on a Riemann surface, 
shown in this figure. The Riemann surface consists of top 
and bottom sheets, which are connected along the segment 
(—00,0] of the real axis (the cut of the Riemann surface into 
sheets is not unique, but this is the standard cut adopted by 
most scientific softwares). Given a point z in the complex 
plane, we can assign to it two points (see P and R in the di- 
agram) on the Riemann surface. Most scientific software will 
automatically assign point P and compute the square root 
of z accordingly. If point R would have been assigned, then 
the square root would differ by exactly a factor —1. When 
computing the square root of a z that is continuously varied, 
one can easily map the position of z on the Riemann surface 
and therefore correct the value computed by the software. As 
examples, we considered three paths for starting at the 
solid circle and ending at the solid square. The square root 
at the end of each path must be corrected by: a = +1 (no 
correction) for path 1 and a = —1 for paths 2 and 3. 



Even for = 0, when inversion symmetry is present, 
the model displays a fairly complex phase diagram as 
function of A. The energy gap closes at A = 0, 27.44, 
116.44, 140.14, 178, 205.44, and the model displays 
QSH (Z2 = -1) phases inside the intervals: (0,27.44), 
(116.44,140.14) and (178,205.44). This can be verified 
directly by computing the parities at the time-reversal 
invariant /c-points. We have verified that our formula for 
the strong Z2 invariant, given in Eq. 44, gives the same 
results. We will not present these calculations here and 
instead we will present in detail the case when the in- 
version symmetry is absent. For this we chose = 2, in 
which case the phase diagram as function of A is qualita- 
tively changed, by the emergence of few metallic phases. 
We will focus only on the lower part of the diagram, 
where direct bands structure calculations show that the 
energy gap closes at A = 10 and that it remains closed 
until A = 13. After that the gap stays opened until 
A = 30, when the gap closes and remains closed when 
A is further increased. Additional phases emerge after 
that, but will not be discussed here. The numerical re- 
sults are shown in Fig. 3, where the two insulating phases 
mentioned above are sampled in 3 points. The calcula- 
tion shows that the insulating phase below A = 10 is 
trivial, while the one between A = 13 and A = 30 is 
topologically non-trivial. 



V. EXTENSION TO CONTINUUM MODELS 

Let us consider a periodic crystal described by a Hamil- 
tonian {II = a lattice vector): 



with (eV units are assumed): 

M = 13.72 cos + 89 cos ^/k^ - + A - 102.72 
Ai = 2.26 sin k^^ A2 = 3.33(sin kx — i sin ky). 

(46) 

The parameter A was allowed to vary. 

The action of the time-reversal operation in is the 
same as described in Eq. 37 and the inversion operation 
is implemented by: 



(47) 



As written above, the model has both time-reversal and 
inversion symmetries. Therefore, we introduce an addi- 
tional term in the Hamiltonian, which breaks the inver- 
sion symmetry: 
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(48) 



iJ = -V^ + t>(r), V{r + R) = V{r) 



(49) 



acting on the Hilbert space H of 2-component spinors 
that are square integrable. The Hamiltonian is assumed 
to commute with the time-reversal operation e^^^^K {K 
= complex conjugation). 

We now consider the Bloch decomposition, given by 
the isometry: 



(50) 



where H is the original Hilbert space and H represents 
the space of square integrable spinors defined over only 
one unit cell. Under this isometry, we have: 



uHu 



(51) 



where Hj^ is given by — + ^(^) but this time defined 
only over one unit cell and with the Bloch boundary con- 
ditions (the prime indicates the derivative): 



^fc(r + i^) = e^'^•^^fc(r) 



ik R 



(52) 
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FIG. 3. Results for the 3- dimensional model of Eqs. 45, 46 and 48, with R — 2. Each panel contains two plots, showing the 
path in the complex plane of det{[/(o,fcy,-7r),(o,fc^,7r)} (left plot), and of det{[/(7r,fcy,-7r),(7r,fcy,7r)} (right plot), both as functions 
of /cy, which is varied from to tt. The values of A and of the resulting a and 'E2D are also shown. The strong invariant 'E^d 
is computed to be +1 for panels (a) and -1 for panels (b). The determinants at /ccc = tt seem to be pinned at +1, which is a 
peculiarity of the model. 



whenever r and r ^ R are on the boundaries of the unit 
cell. The time-reversal operation satisfies: 

OHkO-^=H_k, (53) 

so at this point there is no practical difference between 
the continuum model and the tight-binding models dis- 
cussed in the previous sections. 

The projectors Pk of the Hamiltonians onto the 
states below a given Fermi level are routinely computed 
by the first principle codes. Hence, the monodromies, the 
E2D and EsD invariants can be computed in a straightfor- 
ward fashion. We are currently working on implementing 
the whole construction in our first principle code but, un- 
fortunately, we cannot show any concrete results at this 
time. 



VI. CONCLUSIONS 

Using a monodromy technique, we proposed new for- 
mulations of the Z2 invariants for topological insulators 
with time-reversal symmetry. The formulations are man- 
ifestly gauge independent and we argued that they can 
be effortlessly integrated in the tight-binding as well as 
first principle simulations. Test calculations confirmed 
a full agreement between the new formulations and the 
already established ones. We hope that the expressions 
of the Z2 invariants given in this work will help the sci- 
entists searching for novel non centro-symmetric 3 di- 
mensional topological insulators. We are currently in- 



tegrating the new formulations of the Z2 invariants in 
our first principle electronic structure codes as an au- 
tomated post-processing routines. We hope that other 
electronic structure practitioners will follow our example. 
One other hope of ours is that this gauge independent 
formulations will lead to more effective and transparent 
real space formulations of the Z2 invariants, absolutely 
necessary for understanding the disorder effects in time- 
reversal invariant topological insulators. 

At the end, let us comment about the monodromy 
technique. It definitely helped us avoid complex calcu- 
lations, as the present results followed entirely from the 
group property and the behavior under time-reversal of 
the monodromy. So far, the monodromy technique has 
been applied to inversion symmetric insulators, to fila- 
mentary structures with inversion symmetry supporting 
topological phonon modes, and to time-reversal sym- 
metric insulators. It is very likely that other point sym- 
metries could be handled in a similar fashion, which is a 
future direction that we are currently exploring. 
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